; docformat = 'rst'
;
;+
;
; :Purpose:
;   Run 12-box model using a given set of input parameters
;
; :Inputs:
;
; :Requires:
;
; :Keywords:
;
; :Outputs:
;
; :Example::
;
; :History:
; 	Written by: Matt Rigby, University of Bristol, Aug 15, 2013
;
;-
function a2i_model_run, parameters, measurements, model, lifetime, restart=restart, $
  pollutant=pollutant, startT=startT, endT=endT, in_Restart=in_restart, $
  x=x, state_element=x_element, out_lifetime=out_lifetime, burden=burden, all=all

  compile_opt idl2, hidden

  if n_elements(pollutant) gt 0 then begin
    piMin=pollutant
    piMax=pollutant
  endif else begin
    piMin=0
    piMax=parameters['NPOLLUTANTS']-1
  endelse

  out_lifetime=fltarr(parameters['NPOLLUTANTS']*parameters['TIMESIZE'])
  lifetime_skipped=intarr(parameters['NPOLLUTANTS']*parameters['TIMESIZE'])

  restart=hash()
  burden=hash()

  if keyword_set(all) then begin
    y_out=hash()
  endif else begin
    y_out=dblarr(measurements['NMEASUREMENTS'])
  endelse
  y_skipped=dblarr(measurements['NMEASUREMENTS'])

  for pi=piMin, piMax do begin

    ;DECIDE WHETHER TO RUN POLLUTANT OR NOT
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    do_pollutant=1
    if n_elements(startT) gt 0 then begin
      if startT ge measurements['START', pi] then begin
        startT_model=startT - measurements['START', pi]
        endT_model=endT - measurements['START', pi]
      endif else begin
        ;Check if invariant parameter (e.g. IAV=0 in which case x['TI'] LT 0)
        if x['TI', x_element] lt 0 then begin
          startT_model=0L
          endT_model=endT - measurements['START', pi]
        endif else begin
          do_pollutant=0
          y_skipped[where(measurements['POLLUTANT'] eq pi, count)]=1
          if count eq 0 then stop
          lifetime_skipped[where(lifetime['POLLUTANT'] eq pi, count)]=1
          if count eq 0 then stop
        endelse
      endelse
    endif

    if n_elements(x_element) gt 0 then begin
      if X['ID', x_element] eq 'PO' then begin
        if parameters['OH_A', pi]*exp(parameters['OH_ER', pi]/293.) lt 1.e-16 then begin
          do_pollutant=0
          y_skipped[where(measurements['POLLUTANT'] eq pi, count)]=1
          if count eq 0 then stop
          lifetime_skipped[where(lifetime['POLLUTANT'] eq pi, count)]=1
          if count eq 0 then stop
        endif
      endif
    endif


    ;RUN MODEL
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    if do_pollutant then begin

      if keyword_set(in_restart) then begin
        in_restart_temp=in_restart[pi]
      endif else begin
        in_restart_temp=!null
      endelse

      ;Use OH file (optional)
      OH_input_extension=''
      OH_start_year=-1
      if parameters.haskey('OH_FILE') then begin
        if parameters['OH_FILE'] ne '' then begin
          OH_input_extension=parameters['OH_FILE']
          calendar_year=float(measurements['START', pi])/12. + float(parameters['STARTY'])
        endif
      endif
      ;Use transport file (optional)
      ;IAV NOT YET IMPLEMENTED
      input_extension=''
      if parameters.haskey('TRANSPORT_FILE') then begin
        if parameters['TRANSPORT_FILE'] ne '' then begin
          transport_input_extension=parameters['TRANSPORT_FILE']
          calendar_year=float(measurements['START', pi])/12. + float(parameters['STARTY'])
        endif
      endif

      if parameters['IAV'] then begin
        T_scale=model['T_SCALE', *, measurements['START', pi]:-1]
      endif else begin
        T_scale=model['T_SCALE']
      endelse

      ;Run model
      conc=agage_box_model(model['IC', pi], model['EMISSIONS', pi], $
        arr_OH=[parameters['OH_A', pi],parameters['OH_ER', pi]], mol_mass=parameters['MOL_MASS', pi], $
        lifeTime=model['LIFETIME', pi], T_scale=T_scale, $
        oh_scale=model['OH_SCALE', *, measurements['START', pi]:*], $
        restart=restart_temp, in_restart=in_restart_temp, startT=startT_model, endT=endT_model, $
        overall_lifetime=lifetime_temp, burden=burden_temp, all=all, $
        transport_input_Extension=transport_input_Extension, OH_input_extension=OH_input_extension, $
        calendar_year=calendar_year, IAV=parameters['IAV'])

      restart[pi]=restart_temp
      burden[pi]=burden_temp

      ;Put into output array
      if keyword_set(all) then begin
        y_out[pi]=conc
      endif else begin
        for bi=0, 3 do begin
          y_out[where( (measurements['BOX'] eq bi) and $
            (measurements['POLLUTANT'] eq pi) )]=conc[bi, *]
        endfor
      endelse

      ;OUTPUT INVERSE LIFETIME
      wh=where(lifetime_temp ne 0., count)
      if count eq 0 then stop
      lifetime_temp[wh]=1./lifetime_temp[wh]
      out_lifetime[where(lifetime['POLLUTANT'] eq pi and lifetime['TI'] ge measurements['START', pi])]=lifetime_temp

    endif ;else begin  ;Do pollutant

  endfor  ;POLLUTANT

  if not keyword_set(all) then begin
    for pi=0, parameters['NPOLLUTANTS']-1 do begin
      case strUpCase(parameters['FUNCTION_Y', pi]) of
        'LN(Y)': begin
          wh=where(y_out gt 0. and measurements['POLLUTANT'] eq pi, count)
          if count gt 0 then y_out[wh]=alog(y_out[wh])
        end
        'Y': ;Do nothing
      endcase
    endfor
    
    ;Replace skipped pollutants with reference values
    ;Must come after function_y
    if total(y_skipped) gt 0. then begin
      wh=where(y_skipped gt 0)
      y_out[wh]=measurements['MODEL', wh]
    endif
  endif

  ;Replace skipped pollutants with reference values
  ;Must come after function_y
  if total(y_skipped) gt 0. then begin
    wh=where(lifetime_skipped gt 0, count)
    if count eq 0 then stop
    out_lifetime[wh]=lifetime['MODEL', wh]
  endif
  
  
  return, y_out

end